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Abstract 

A phylogenetic tree and median-joining network based on cytochrome b sequence data revealed clades consistent with 
morphological differences and geographical distribution of Clarias batrachus (Linnaeus, 1758) in Southeast Asia. AMO- 
VA analysis for variation was significant among populations (/ > <0.05) and was in agreement with morphological differ¬ 
ences. Pairwise differences were significant between Java and Brunei/Bomeo, Brunei/Bomeo and west Malaysia, and 
Java and west Malaysia samples (P < 0.05). Closest relationships were found between samples from Brunei/Borneo and 
Java, and between west Malaysia and Laos-Sumatra. Nine haplotypes were unique to geographical regions. The Java spe¬ 
cies had high haplotype (1.000 ± 0.126) but low nucleotide (0.017) diversities, suggesting a population bottleneck fol¬ 
lowed by expansion. However SSD and Hri {P= 0.5) did not support demographic expansion. Instead, purifying selection 
where mutations occur and accumulate at silent sites is a more acceptable explanation. 
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Introduction 

In Southeast Asia, 19 species of the genus Clarias are currently recognized (Ng et al., 2011) including 13 species in 
Borneo (Kottelat, 2013). In spite of the relatively recent colonisation of Southeast Asia from central Asia about 15 
MYA (Agnese & Teugels, 2005), Clarias shows considerable diversity in morphology and habitat. 

Clarias batrachus is a common species found in South and Southeast Asia and is widely utilised in aquaculture 
and the aquarium fish trade, and has been the subject of numerous studies (see Begum, 2004; Hora, 1936; Kottelat, 
2001; Ng & Kottelat, 2008; Khedkar et al. 2010). The species has been introduced to many parts of the world and 
has the potential to threaten non-native environments (Lever, 1996). 

The identity of C. batrachus has recently been studied by Ng & Kottelat (2008). They redescribed C. 
batrachus and designated a neotype from Java. Ng & Kottelat (2008) designated the non-Javanese C. batrachus as: 
Clarias aff. batrachus ‘Indochina’ from the Mekong River drainage, and Clarias aff. batrachus ‘Sundaland’ from 
west Malaysia and Sarawak. They noted morphological differences in these forms in the shape of the supraoccipital 
process and width of the frontal fontanelles. The morphological difference between C. aff batrachus ‘Sundaland’ 
and C. aff. Indochina was in the shape of the supraoccipital process. The latter was more triangular. The shape of 
supraoccipital process also is a useful diagnostic character in Asian Clarias (Teugels et al. 1999). The frontal 
fontanelle of C. batrachus from the Malay Peninsula and Sarawak was wider and referred to as ‘sole shaped’ by 
Teugels (1986) than that on specimens from Java and mainland Southeast Asia, referred to as ‘knife-shaped’ by 
Teugels (1986). The morphological differences seemed to correlate with geographical areas of Southeast Asia. 

Several molecular genetics studies had been conducted on Asian Clarias (Agnese & Teugels, 2005; Pouyaud et 
al. 2009), Asian C. batrachus (Islam et al. 2007; Khedkar et al. 2010) and Southeast Asian C. batrachus (Mohindra 
et al. 2007). However, no molecular study has investigated variation in sequence data from the cytochrome b gene 
in C. batrachus from Southeast Asia. The objective of the present study is to determine whether genetic data from 
the cytochrome b gene is in agreement with morphological variation in C. batrachus from Southeast Asia. 
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FIGURE 1. Dorsal views of heads of: A, Clarias batrachus from Brunei/Borneo, UBD 0412, 202.6 mm SL; B, C. batrachus 
from Java, UBD 0814, 221.8 mm SL. Frontal fontanelles ‘sole-shaped’ in A and ‘knife-shaped’ in B; less triangular shape of 
the supraoccipital process in both A and B. 
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FIGURE 2. Maximum likelihood (GTR) tree of Clarias batrachus cytochrome b gene. Numbers on the nodes correspond to 
branch length/ML bootstrap probability except * represents NJ bootstrap value/ branch length. Only values >50% are shown. 
Brunei/Borneo (BN1-15), Java (JAV1, 2, 4, 5, 6), west Malaysia (JF.MY1-5), Laos (AB.LAOI), and Sumatra (AF.SUM1). 
Morphological characters are labeled on each Clade. 


184 ■ Zootaxa 3962 (1) © 2015 Magnolia Press 


LEE & SULAIMAN 



























Materials and Methods 


Fish samples and morphology 

Specimens of C. batrachus (Fig.1) were collected from the Tutong District in Brunei, Borneo (n=15) and east Java, 
Indonesia (n=5). The identity of C. batrachus from Java was confirmed by the recent designation of the type 
locality by Ng & Kottelat (2008). Fish sampling in Brunei was done using gill nets and fishing rods. Samples from 
Indonesia were obtained from a local fish market at Bekasi, Java. The fishes were stored at -20°C until subjected to 
DNA extraction. Material examined in this study is deposited at the Universiti Brunei Darussalam, Brunei (UBD). 


DNA extraction and sequencing 

Total genomic DNA was extracted from 30-100 mg tissue using DNeasyTM tissue kits (Qiagen GmBH, 
Germany). Amplifications were carried out in 100 pi containing lx reaction Buffer, 3.5mM MgCl 2 , 0.2mM dNTP, 
1 pM of forward and reverse primers, 2.5 units of Tag DNA polymerase (Invitrogen, San Diego, CA) and about 
50-100 ng of genomic DNA. The 435 bp partial sequences of mitochondrial cytochrome b gene were amplified 
with primers NEW-FOR 5’ AGC CTA CGA AAA ACC CAC CC 3’ (Chang et al. 1994) and H15149 5’ AAA 
CTG CAG CCC CTC AGA ATG ATA TTT GTC CTC A 3’ (Kocher et al. 1989). 

Thermal cyclings were performed using GeneAmp PCR system 9700 (Applied Biosystem) with an initial 
denaturation of 94°C for 5 minutes, followed by 30 cycles of 94°C at 30 seconds, annealing temperature of 55°C 
for 30 seconds, 72°C for 60 seconds and the final extension of 5 minutes at 72°C. PCR products were purified 
using a Qiagen purification kit according to manufacturer’s instructions and subsequently sequenced at First Base 
(www.base-asia.com). 


Data analysis 

All 27 sequences were edited using Clustal X version 2.0 (Larkin et al. 2007). Missing nucleotide data were 
excluded on a pair-wise basis, ffaplotypes (H). haplotype diversity ( Hd) (Nei, 1987) and nucleotide diversity (/;) 
(Nei, 1987) were calculated for each geographic population using DnaSP 5 (Librado & Rozas, 2009). Nucleotide 
composition, parsimony informative and variable sites were estimated using MEGA 6.06 (Tamura et al. 2013). 

A maximum-likelihood (ML) tree based on General Time Reversible (GTR) model (Nei & Kumar, 2000) and 
a neighbour-joining (NJ) tree based the Tamura 3-parameter (Tamura, 1992) were generated using MEGA 6.06 
(Tamura et al. 2013). Trees were unrooted to specify only relationships among C. batrachus. The robustness of 
statistical support for branches on the ML and NJ trees was determined by 1000 bootstrap replicates (Felsenstein, 
1985). The branch length measured the number of substitutions per site. The evolutionary relationship among 
cytochrome b haplotypes was resolved for C. batrachus with unrooted networks constructed with the program 
Network 4.6 (http://www.fluxus-engineering.com) using a median-joining algorithm and default settings (Bandelt 
et al. 1999). The partial cytochrome b gene sequences from west Malaysia (GenBank Accession Nos. JF280847, 
JF280849, JF280854, JF280855 and JF280859), Laos (GenBank Accession No. AB8222528), and Sumatra 
(GenBank Accession No. AF 235932) were included in the trees and the MJ network. 

An analysis of molecular variance (AMOVA; Excoffier et al., 1992) was conducted in Arlequin 3.0 (Excoffier 
et al., 2005) to test the spatial genetic variation of C. batrachus species within the five areas. Additionally, the 
fixation index F SJ value was calculated from haplotype frequencies, to determine the level of genetic differentiation 
(Weir & Cockerham, 1984). Tajima’s D (Tajima, 1989a; 1989b) and Fu’s F (Fu, 1997) statistics were also 
conducted through Arlequin 3.0 (Excoffier et al., 2005) as well to understand the degree of population expansion 
and to test the deviation from equilibrium expectations. The goodness of fit test between the observed and expected 
distributions was estimated based on the JJarpending’s raggedness index (Hri) and the sum of squared deviations 
(,SSD ) for the expected stepwise expansion model (Harpending, 1994). The sudden expansion model 6 and z 
statistics parameters were used to measure changes in population size (whereas 6 was being expressed in a 
mutational and general scale) and the age of populations in mutational time was calculated to test whether the 
observed mismatch distributions in each population better fit the model of constant population size or the one 
model with population size changes (sudden expansion model) (Rogers, 1995). 
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Results 


Molecular characteristics 

Alignment of 27 sequences of the cytochrome b gene resulted in 395 sites, of which 74 were variable (18.7%), and 
61 were parsimony informative (15.4%). Nine haplotyes were identified. The base composition was 32.4% (A), 
29% (G), 21.8% (C) and 16.8% (T). The Java specimens were the most diverse (Haplotype diversity, Hd was 1.000 
± 0.126), and nucleotide diveristy (ir) was 0.01722 (Table 1). The Brunei/Borneo and west Malaysia populations 
had 1 haplotype each and 0.00 for haplotype and nucleotide diversities. Tajima’s D and Fu’s F for C. batrachus 
from Java were -0.525 and -0.145, respectively, and not statistically significant. 


TABLE 1. Sampling regions, number of samples (N), number of haplotypes ( H ), haplotype diversity (Hd), nucleotide diversity 
(p), Tajima’s D and Fu’s Fs for cytochrome b gene sequences for each population of Clarias batrachus. 


Sampling region 

N 

H 

Hd 

p 

Tajima's D 

Fu's Fs 

East Java/Bekasi 

5 

5 

1.00±0.126 

0.01722 

0.08677 

-0.875 

Malaysia 

5 

1 

0.000 

0.000 

0.000 

0.000 

Brunei 

15 

1 

0.000 

0.000 

0.000 

0.000 

Laos PDR 

1 

1 

- 

- 

- 

- 

Sumatra 

1 

1 

- 

- 

- 

- 

Total 

27 

9 






Patterns of population structure 

The phylogenetic analyses revealed four geographically exclusive clades (Fig. 2). The topology of ML and NJ 
were similar. The bootstrap probabilities and branch lengths from ML and NJ were in agreement except where it 
was marked as * on the ML tree. Clade A contained 100% of C. batrachus from Brunei/Borneo and was supported 
by 99% (NJ)/ 52% (ML). Clade B contained 100% of C. batrachus from Java and supported by 100% bootstrap. 
Clade C contained specimens from Laos and Sumatra, and was supported by 100% bootstrap. Clade D contained 
100% of the C. batrachus from west Malaysia and was supported by 100% bootstrap. Clades A-B and C-D were 
resolved as sister taxa with 100% bootstrap. Evolutionary relationships using a median-joining network were 
further analysed between haplotypes of C. batrachus. The haplotype networks (Fig. 3) supported the topology of 
the phylogenetic tree in Figure 2. It formed four haplotype groups (A, B, C and D). Flaplotype groups A (Brunei/ 
Borneo) and B (Java) were separated by 13 mutational steps. Specimens from Java showed five haplotypes and 
could represent various populations. A and D (west Malaysia) haplotype groups were separated by 38 mutational 
steps. Haplotype groups C split into (i) Sumatra and (ii) Laos, and each was separated from west Malaysia by 6 
mutational steps. 

The pairwise 7%^ values ranged from 0-1.0 and agreed with the phylogenetic tree. It was significant only 
between samples from Java and Brunei/Borneo, C. batrachus from Java and west Malaysia, and C. batrachus from 
Brunei/Borneo and west Malaysia (Table 2). Genetic differentiation was high, indicating low gene flow and 
connecivity among populations, except pairwise differences between Java-Sumatra and Java-Laos were 0.0. 

AMOVA analysis showed that genetic variation observed in C. batrachus was significant (7 > <0.05) among 
populations (Table 3). The tests for SSI) and Hri were performed on C. batrachus from Java only as mismatch 
variances for other samples were too small. It had a multimodal and rough mismatch distribution (Fig. 4), 
consistent with old stationary populations. Flri was higher (0.2) and SSD was lower (0.06), but both tests were 
insignificant (/*=().5), suggesting a poor fit to demographic expansion. 
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FIGURE 3. Median-joining network derived from Clarias batrachus cytochrome b sequences. Each circle represents a 
haplotype and is scaled to the number of individuals with that haplotype. Small dark spots represent mutational steps between 
haplotypes. 
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FIGURE 4. Mismatch distribution based on cytochrome b gene sequences of C. batrachus from Java. Expected distribution 
(Exp) and Observed distribution (Obs) under the model of population expansion. 


TABLE 2. Population pairwise F SJ values for cytochrome b data of C. batrachus are given below, and F ST p values are given 
above the diagonal. Java (JAV), Brunei (BN), west Malaysia (MY), Sumatra (SUM), and Laos (LAO). Significant p values are 
indicated by an asterisk. 



JAV 

BN 

MY 

SUM 

LAO 

JAV 

0.0000 

<0.05* 

<0.05* 

0.999 

0.999 

BN 

0.71963 

0.0000 

<0.05* 

0.999 

0.999 

MY 

0.50000 

1.0000 

0.0000 

0.999 

0.999 

SUM 

0.00000 

1.0000 

1.0000 

0.0000 

0.999 

LAO 

0.00000 

1.0000 

1.0000 

1.00000 

0.0000 
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TABLE 3. AMOVA analysis for C. batrachus in Southeast Asia. 


Source of variation 

d.f. 

Sum of squares 

Variance components 

Percentage of variation 

Among populations 

4 

6.741 

0.381Va (p<0.05) 

80.73 

Within populations 

22 

2.000 

0.091 Vb 

19.27 

Total 

26 

8.741 

0.472 



Fixation Index F S1 : 0807 


Discussion 

In this study, base compositions were similar to those found in Asian Clarias (Pouyaud et al. 2009) and in 
agreement with the pattern described for numerous fishes (Cantatore et al. 1994). 

The Brunei/Borneo-west Malaysia pair had the highest number of distinct mutations, followed by the Brunei/ 
Borneo-Java pair. This was supported by significant pairwise F ST values between these populations (P<().()5). The 
high genetic differentiation was expected as freshwater fishes are confined to specific basins (Ward et al. 1994). 

The genetic variation in this study was consistent with morphological variation in the supraoccipital process 
(less triangular in Java), frontal fontanelle (‘knife-shaped’ in Java and ‘sole shaped' in Brunei/Borneo). 

The haplotypes from the cytochrome b gene data identified distinct groups of C. batrachus from Java, Brunei/ 
Borneo, west Malaysia and Lao-Sumatra. The single haplotype for Brunei/Borneo was the result of two mutations 
(Fig. 3). The Laos and Sumatra specimens split from the Malaysian population to form sister branches. Ffowever, 
C. batrachus was represented from Laos and Sumatra by one specimen each, and results obtained in this study 
might not be representative of the populations. More samples will be required to confirm this relationship. 

The current study showed a lack of genetic differentiation between the Java-Laos pair and the Java-Sumatra 
pair. The latter may be caused by the close proximity of the North and East Sunda paleo-drainages (Rainboth, 
1991). The results for lack of genetic differentiation between Java and Laos pair was consistent with a “radiation by 
paleo-drainage” model of diversification observed for halfbeaks across Southeast Asia (De Bruyn et al. 2013). 

High haplotype and low nucleotide diversities of C. batrachus from Java suggested a rapid population growth 
from an ancestral population with small Ne after a bottleneck event (Avise, 2000). The results supported the 
hypothesis of a pre-Pleistocene origin in Asian Clarias diversity, with colonization of Southeast Asia from Central 
Asia about 15 MYA (Agnese & Teugel, 2005) and appeared on the Sunda Shelf (Sumatra, Borneo and Java) about 
20 MYA (Pouyaud et al. 2009). The negative D <0.1 may be interpreted as due to increasing population size or as 
purifying selection at a locus (Tajima 1989a). Negative F means a recent population expansion or purifying 
selection. However, population expansion of the Java population was not supported by a mismatch distribution, F[ri 
and SSD (P >0.1). Purifying selection represents the likely alternative, as shown by the negative D and F values. 
Purifying selection occurs when mutations accumulate at silent sites but with little heterozygosity (Tajima, 1989a). 

Our study revealed significant molecular diversification within C. batrachus. These results are consistent with 
morphological disparity and reproductive isolation of the Southeast Asian C. batrachus populations. Greater taxon 
and character sampling throughout the range (especially Indo-China and Sumatra) of C. batrachus is necessary to 
better determine the variation and species boundaries of C. batrachus. 
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